Abstract. We present the results of three-dimensional 
hydrodynamical simulations of the subsonic thermonu- 
clear burning phase in type la supernovae. The burning 
front model contains no adjustable parameters so that 
variations of the explosion outcome can be linked directly 
to changes in the initial conditions. In particular, we in- 
vestigate the influence of the initial flame geometry on the 
explosion energy and find that it appears to be weaker 
than in 2D. Most importantly, our models predict global 
properties such as the produced nickel masses and ejecta 
. velocities within their observed ranges without anv fine 
^ _ tuning. 
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1. Introduction 



In a series of papers ( Reinecke et al. 1999b| ,a|, |2QQ2|) we 
described a new numerical tool to simulate thermonuclear 
explosions of white dwarfs in two and three spatial di- 
mensions. Our aim was to construct models of type la 
supernova (SN la) explosions that are as free from non- 
physical parameters as currently feasible. Thermonuclear 
burning, in particular, is represented by a subsonic turbu- 
lent flame whose local velocity is derived from a subgrid- 
scale model for unresolved turbulent fluctuations. Solved 
in combination with the compressible Euler equations, this 
model contains no free parameters that could be adjusted 
in order to fit SN la observations. Consequently, the ini- 
tial white dwarf model (composition, central density, and 
velocity structure), as well as assumptions about the lo- 
cation, size and shape of the flame surface as it first forms 
fully determine the simulation results. Here, we concen- 
trate on variations of the latte r. 

In Reinecke et al. (1999b|), we confirmed the earlier 
result of |Niemeyer et al. (1996| ) that, at least in 2D, the 
explosion energy and amount of 56 Ni produced in the ex- 
plosion (which determines the brightness of an SN la) are 
sensitive to the ignition conditions. To be more precise, 
a more complicated topology of the initial nuclear flame 
seems to lead to higher Ni-production and, consequently, 
more powerful explosions. One might even speculate that 
the randomness of the ignition process could be the rea- 
son fortheobseryed^ of normal SNe 



la ( [Hillebrandt fc Niemeyer 2000D . 

This article continues the presentation of numerical 
simulations of SNe la by testing the effect of different ini- 
tial conditions on the simulation outcome in three dimen- 
sions. In this context, the simultaneous runaway at several 
different spots in the central region of the progenitor star 
is of particul ar interest. A plausible ignition scenario was 
suggested by Garcia-Senz fc Woosley (1995 ). 

All simula tions were performed u sing th e algorithms 
presented by [Reinecke et al. (2QQ2| , |l999b|) . The initial 
model for the white dwarf, assumed to be at the Chan- 
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drasekhar mass, as well as the grid geometry and symme- 
try assumptions, are identical to the setup described in 
Sects. 3 and 4 of |Reinecke et al. (200 j) . 

Sect. U of this paper presents the initial front geom- 
etry and its temporal evolution for two multi-point igni- 
tion calculations. Various aspects of all three-dimensional 
calculations are then compared in Sect. [| As far as pos- 
sible, this comparison is then extended to other existing 
results of SN la simulations. In the case of parameterized 
one-dimensional calculations such comparisons are diffi- 
cult, since most of the published data have no analogy in 
three dimensions. As an alternative, we suggest to track 
and compare the amount of burned material as function 
of the density at which the burning took place. 

Generally speaking, we confirm the earlier 2D results: 
models with more, or better resolved, ignition spots tend 
to produce more radioactive Ni also in 3D, although the 
effect is somewhat smaller. But together with the gain of 
Ni by going from 2D to 3D our models predict amounts 
that are in g ood agreement with t hose inferred from the 
observations ( Contardo et al. 2000 ). This will also be dis- 
cussed in Sect. |3|, and our conclusions follow in Sect. O. 



2. Multi-point ignition scenarios 

Two different calculations were performed to investigate 
the off-center ignition model in three dimensions. The sim- 
ulation b5_3d_256 was carried out on a grid of 256 3 cells 
with a central resolution of 10 6 cm and contained five bub- 
bles with a radius of 3 • 10 6 cm, which were distributed 
randomly in the simulated octant within 1.6 • 10 7 cm of 
the star's center. As an additional constraint, the bubbles 
were required not to overlap significantly with the other 
bubbles or their mirror images in the other octants. This 
maximizes the initial flame surface. The exact algorithm 
for the positioning of the bubbles was the following: First, 
the centers of the bubbles were chosen randomly within 
the radius of 1.6 • 10 7 cm. A realization was accepted when 
the distance between any two bubble centers (including 
the mirror images across the coordinate planes) was larger 
than 1.3 bubble radii. A particular realization of these ini- 
tial conditions and its temporal evolution is shown in Fig. 
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Fig. 1. Snapshots of the flame front for scenario b5_3d. 
The fast merging between the leading and trailing bubbles 
and the rising of the entire burning region is clearly visible. 
One ring on the coordinate axes corresponds to 10 7 cm. 



Part. No. 


1 
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3 


4 


5 


x-Pos. 


128 


35 


42 


85 


63 


y-Pos. 


76 


60 


129 


27 


85 


z-Pos. 


27 


30 


30 


41 


83 



Table 1. Approximate position of the bubble centers in 
the model b5_3d_256. All lengths are given in kilometers. 



1. The locations of the bubble centers are listed in Table 
1. 

In a way very similar to the evolutio n of earlier two- 
dimensional simulations (cf. Fig. 6 of Reinecke et al. 
1999a0 the flame kernels closer to the center are elon- 
gated very rapidly and connect to the outermost bubbles 
within 0.15 s. The whole burning region disconnects from 
the coordinate planes and starts to float slowly towards 
the stellar surface. 

In an attempt to reduce the initially burned mass 
as much as possible without sacrificing too much flame 
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Fig. 2. Snapshots of the flame front for the highly resolved 
scenario b9_3d. One ring on the coordinate axes corre- 
sponds to 10 7 cm. 



surface, the very highly resolved model b9_3d_512 was 
constructed. It contains nine randomly distributed, non- 
overlapping bubbles with a radius of 2-10 6 cm within 
1.6 • 10 7 cm of the white dwarf's center. To properly rep- 
resent these very small bubbles, the cell size was reduced 
to A = 5 • 10 5 cm, so that a total grid size of 512 3 cells 
was required. 

Starting out with very small flame kernels is desir- 
able, because the initial hydrostatic equilibrium is bet- 
ter preserved if only a little mass is burned instanta- 
neously. Furthermore the floating bubbles are expected 
to be even smaller in reality than in the presente d model 
(r J$ 5 • 10 5 cm, see |Garcia- Senz fc Woosley 1995|) . 

The center locations of our particular realization are 
given in Table ||. Snapshots of the front evolution (Fig. 
^ exhibit features very similar to those observed in the 
lower-resolution 3D models. Only in the last plot the for- 
mation of additional small-scale structures becomes evi- 
dent. 
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Model b5_3d_256: composition evolution 
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Fig. 3. Time evolution of the chemical composition in the models b5_3d_256 and b9_3d_512. 
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x-Pos. 


47.7 


41.4 


27.5 


42.2 




y-Pos. 


45.5 


133.1 


24.1 


85.0 
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34.0 


39.3 


21.7 
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Energy comparison, 3D simulations 



Table 2. Position of the bubble centers in the model 
b9_3d_512. All lengths are given in kilometers. 



The direct comparison of the nuclear evolution with 
model b5_3d_256 in Fig. || nevertheless reveals differ- 
ences in the early explosion stages: the high-resolution 
simulation produces slightly more nickel and considerably 
more a-particles during the first half second. This phe- 
nomenon is most likely explained by the discrepancy in 
the ratio between the initial flame surface and the vol- 
ume of burned material in the two models. The higher 
burned mass in the five-bubble model initiates an early 
bulk expansion of the star and therefore causes a rapid 
drop of the central density. Since its flame surface is quite 
small compared to the burned volume, only relatively little 
mass can be burned at high densities. In the nine-bubble 
model, on the other hand, the star expands more slowly at 
first because less material is burned instantaneously, and 
the mass fraction of a-particles in the reaction products 
is consequently rather high. Since the energy buffered in 
those a-particles is not immediately used to drive the ex- 
pansion, the flame can consume more material at higher 
densities and has more time to increase its surface as a 
result of hydrodynamical instabilities. 

3. Comparison of all 3D simulations 

3.1. Energy release 

In addition to the two simulations presented above, the 
following comparison also includes the model c3_3d_256, 
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Fig. 4. Energy evolution of the three-dimensional explo- 
sion models of Fig. 1 (dashed) and Fig. 2 (dashed-dotted). 
For comparison w e also show the centr ally ignited ("three 
fingers") model of Reinecke et al. (2002| ) (solid line) . While 
the centrally ignited and the five-bubble models are re- 
markably similar, the high-resolution nine-bubble simula- 
tion reaches higher energies despite relatively slow initial 
burning. 



which was described in Reinecke et al. (2002| ). It must 
be emphasized that the comparison between the normal 
and high-resolution calculations can be qualitative at best, 
because the initial conditions as well as the resolution of 
model b9_3d_256 are different from the other two, and it 
was not possible to perform a 3D resolution study like the 
two-dimensional one discussed in Reinecke et al. (200^ ) in 
order to investigate the influence of the cell size on the 
simulation outcome. However, if one assumes that numer- 
ical convergence is reached as long as the flame stays in 
the uniform grid region and the cell size is at most 10 6 cm 
(which seems very plausible according to the 2D study), at 
least the evolution during the initial ^0.8 seconds should 
be practically independent of the employed resolution. For 
example, if the model b5_3d would be re-calculated at in- 
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creased resolution, the authors would expect identical be- 
haviour to the original simulation up to 0.8 seconds after 
ignition, and probably slightly increased energy generation 
afterwards. 

During the first 0.5 seconds, the three models are 
nearly indistinguishable as far as the total energy is con- 
cerned (see Fig. ||), which at first glance appears somewhat 
surprising, given the quite different initial conditions. A 
closer look at the energy generation rate actually reveals 
noticeable differences in the intensity of thermonuclear 
burning for the three simulations, but since the total flame 
surface is initially very small, these differences have no vis- 
ible impact on the integrated curve in the early stages. 

However, after about 0.5 seconds, when fast energy 
generation sets in, the nine-bubble model burns more vig- 
orously due to its larger surface and therefore reaches a 
higher final energy level. Fig. || a lso shows that the ce n- 
trally ignited model c3_3d_256 of Reinecke et al. (2002|) is 
almost identical to the off-center model b5_3d_256 with 
regard to the explosion energetics. 

Obviously, the scatter in the final energies due to differ- 
ent initial conditions appears to be subst antially smaller 



than for simulations in two dimensions (see |Reinecke et al. 
1999b|. This is, in part, good news, since it demonstrates 
the robustness of the explosion mechanism. On the other 
hand, it may lower the chances to explain all "normal" 
type la supernovae by random variations of ignition con- 
ditions alone. Whether this result is a consequence of the 
particular choice of initi al conditions or the se lf-regulation 
mechanism described in Reinecke et al. (2002 ) works more 
efficiently in 3D needs to be studied further. 

3.2. Distribution of the ejecta in real and velocity space 

It is deduced from observations that the angular distri- 
bution of the ejecta is fairly homogeneous in SN la. This 
applies to the material velocities as well as the chemical 
composition. Concerning the ejecta speeds, the performed 
simulations are in good agreement with observations, since 
the subsonic nature of the combustion process allows for 
fast equilibration of pressure fluctuations and therefore 
generates a practically spherical expansion. However, the 
final distribution of the reaction products consists of a few 
large "lumps" of processed material separated by areas of 
unburned carbon and oxygen. Similar resu lts were also ob- 
tained by Niemeyer fc Hillebrandt (1995| ). This apparent 
inconsistency with observations is possibly explained by 
the still limited numerical resolution, which allows only 
a few initial features like bumps in the front or burn- 
ing bubbles to be prescribed. These in turn will dominate 
the large-scale front geometry at late times. In reality the 
number of bubbles is likely to be much larger, so that the 
initial conditions could be isotropic in a statistical sense, 
which would result in a much more uniform ejecta compo- 
sition. The discussed constraints for the initial conditions 
in the models b5 3d and b9 3d also caused the coordi- 



nate planes to be practically free of burned material. As 
a consequence, the most likely macroscopic flow pattern 
that will develop consists of an updraft near the center (?) 
of the octant and downflow at the edges of the computa- 
tional domain. 

The time evolution of line intensities in SN la spectra 
and also the line widths due to the Doppler broadening 
can be used to determine the radial position of different 
elements in the supernova and their velocity. A rough es- 
timate for the velocity distribution of the ejecta for the 
models c3_3d_256 and b5_3d_256 is shown in Fig. |. It 
must be stressed that these graphs do not represent the fi- 
nal ejecta speeds, since 1.5 seconds after ignition — at the 
end of the simulation — the gravitational binding energy 
is still rather high, and leaving the potential well will cause 
some slowdown, especially in the inner regions. Neverthe- 
less the maximum speed and the relative distribution of 
the elements in velocity space allow for some conclusions 
concerning the predictive power of our simulations. 

The ejected material reaches speeds of up to 
12 000 km/s for both initial conditions which is well within 
the observed range of 10 000 to 15 000 km/s of typical 
SNe la. However, the composition of the ejecta in veloc- 
ity space is quite different from results obtained in one- 
dimensional simulations, insofar as both fuel and ashes 
are present at almost all velocities. Such a situation is, by 
definition, impossible in centrally ignited one-dimensional 
model, because everything inside the radius of the flame 
must have been burned and the material further out can- 
not have been processed. Multidimensional calculations, 
in contrast, allow simultaneous burning at many differ- 
ent radii, and fuel and ashes can be stirred by large-scale 
turnover motions. 

Whether the velocity profiles of the current simula- 
tions are compatible with observations has not been thor- 
oughly investigated so far. The large amount of carbon and 
oxygen at low velocities which is typical for the floating- 
bubble calculations (e.g. in the right hand plot of Fig. |H|) 
has not been detected yet in supernova spectra and, there- 
fore, might not be present in real SN la. Whether this fea- 
ture is typical for offcenter ignition scenarios or caused by 
our prescriptions for the initial flame geometry (see also 
Sect. ^2| ) must be investigated in further calculations. On 
the other hand, it is currently conjectured that a certain 
amount of unburned material could remain undetected by 
current spectroscopy near the center of the remnant (P. 
Mazzali, personal communication), so that the off-center 
ignition scenario cannot be firmly ruled out. 

3.3. The final chemical composition 

Our simulations only employed a minimal nuclear reaction 
network which was sufficient for a good approximation 
of the thermonuclear energy release. For that reason, the 
produced abundances cannot be taken literally, but must 



6 



M. Reinecke et al.: Three-dimensional simulations of type la supernovae 




v r [cm/s] v r [cm/s] 



Fig. 5. Distribution of chemical species in radial velocity space for models c3_3d_256 and b5_3d_256 after 1.5 
seconds. While the elements are more or less uniformly distributed in the centrally ignited model, the floating-bubble 
model still contains large amounts of unprocessed material near the center. 



be interpreted carefully when trying to derive light curve 
shapes or spectra from our results. 

In our simulations, the abundance of 56 Ni actually rep- 
resents the sum of all Fe-group elements synthesized dur- 
ing the explosion. Due to the very similar nuclear binding 
energies of all those nuclei, this simplification is justified as 
far as the explosion energy is concerned. To calculate the 
light curve, however, 56 Ni and the rest of the Fe-group 
elements play significantly different roles. While the ra- 
dioactive decay of 56 Ni supplies the energy for the light 
emitted by the supernova ejecta and thereby determine 
the bolometric luminosity, the presence of the other nuclei 
enhances the opacity of the ejecta and causes a broadening 
of the light curve. 

It is expected that "real" 56 Ni accounts for w 70 — 
80% of the Fe-group elements ( Thielemann et al. 1986| ), 
which means that the corrected ^Ni mass is stil l in rather 
good a greem ent with the amounts postulated by Contardo| 
et al. ( |2000|) . However, since the relative mass fraction of 
56 Ni and the total Fe-group elements may also vary in 
space, reliable answers can only be given once detailed 
post-processing using a large reaction network has been 
performed on the existing data. This work is currently in 
progress. 



4. Comparison to previous simulations 

It has been shown above that, as far as comparisons are 
presently possible, our numerical models produce results 
which agree fairly well with real SN la explosions. Unfor- 
tunately, the only parameters we can directly and quan- 
titatively compare to the observations at present are the 
energy release and - to some extent - the nickel masses; all 
other data derived from the simulations (like the scatter 
in the final energies and the amount of intermediate-mass 
elements) are of qualitative nature only. Although they 
agree roughly with expectations, they may not provide a 



criterion for the validity of the model. In order to allow 
for a more detailed comparison with the observations the 
code has to be extended. Apart from post-processing the 
ejecta to produce detailed nuclear abundances, it might 
also be necessary to model the complex physical processes 
taking place in the later SN stages, up to several weeks 
after explosion. 

As an intermediate step, it might therefore be use- 
ful to perform a detailed comparison of our results with 
one of the successf ul phenomenological S N la models 
like Nomoto's W7 (|Nomoto et al. 1984| ). These one- 
dimensional models contain complex nuclear reaction net- 
works and have been used to generate synthetic spectra 
and light curves. A good candidate for such a comparative 
study would be the amount of burned mass as a function of 
the density at which the burning took place. If such graphs 
were similar for our simulations and the phenomenologi- 
cal models, this would argue for the total chemical com- 
position of both models also being fairly equal, which in 
turn would be a strong hint that light curves and spectra 
of both simulations might not be too different. In other 
words, one would check whether characteristic properties 
of the multidimensional numerical experiments agree with 
the results of the phenomenological models. If so, the pre- 
dictions made by the phenomenological simulations could 
be believed to be valid also for our calculations. 

Care must be taken that only quantities are chosen 
for comparison which have a physical meaning in both 
models. For example, the time dependence of the burn- 
ing speed, which is often used to characterize phenomeno- 
logical models and is documented in many publications, 
cannot be used for this purpose, since there exists no sin- 
gle burning speed in multidimensional simulations at any 
given time. In this concrete situation, the time-dependent 
energy generation rate would be the proper choice for a 
comparison, since it is well defined in both one- and multi- 
dimensional scenarios. It would therefore be convenient if 
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the analysis of phenomenological models contained global 
quantities like energy generation rates or density depen- 
dence of burned mass, which can be easily compared to 
multidimensional simulations. 

Khokhlov (2000) describes an alternative approach to 
model SN la explosions and presents results of three- 
dimensional calculations. The numerical models employed 
in his code differ quite strongly in various aspects from the 
schemes in our code. In particular, the simulation is car- 
ried out on an adaptively refined grid, the flame is modeled 
by a reaction-diffusion method, and its propagation veloc- 
ity is determined by the asymptotic rise speed of Rayleigh- 
Taylor bubbles instead of the turbulent velocity fluctua- 
tions. 

The initial setup described in his paper is centrally 
ignited, and no perturbation is applied to the spherical 
flame surface. As soon as nonlinear instabilities have de- 
veloped, however, the explosion progresses in a way which 
looks very similar to the 3D simulations discussed in Sect. 
|2|. This similarity applies to the energy production rate as 
well as the geometrical features developed by the flame. 

It is an encouraging fact that two independent and 
quite different attempts to simulate a SN la without in- 
troducing artificial burning velocities and chemical mixing 
yield fairly comparable results. This shows that the ex- 
plosion mechanism analysed here in detail, namely a pure 
deflagration front in a C+O Chandrasekhar-mass white 
dwarf, is indeed robust. 

5. Discussion and outlook 

In this paper we have presented the results of three- 
dimensional numerical simulations of thermonuclear defla- 
gration fronts in Chandrasekhar mass white dwarfs com- 
posed of equal amounts of carbon and oxygen. We could 
show that independent of the details of the ignition pro- 
cess (which is still far from being well understood) the 
white dwarf is always disrupted by the release of nuclear 
energy. As far as we could check at present, the properties 
of our models are in good agreement with observations of 
typical type la supernovae. In particular the explosion en- 
ergy and the average chemical composition of the eject a 
seem to fit the observations (see also Table 1). This success 
of the models was obtained without introducing any non- 
physical parameters, but just on the basis of a physical and 
numerical model of subsonic turbulent burning fronts. We 
also stress that our models give clear evidence that the 
often postulated deflagration-detonation transition is not 
needed to produce sufficiently powerful explosions. 

There are certainly several desirable additions and im- 
provements. The most crucial question still seems to be 
the ignition process, although in 3D the dependence of the 
final outcome is weaker than in our previous 2D models. 
In principle one could try to simulate the ignition phase 
with the numerical models we have developed, but be- 
cause of the much longer time scales this requires huge 



model name 


m Mg [M Q ] 


m N i [M Q ] 


E nuc [10 50 erg] 


c3_3d_256 


0.177 


0.526 


9.76 


b5_3d_256 


0.180 


0.506 


9.47 


b9_3d_512 


0.190 


0.616 


11.26 



Table 3. Overview over element production and energy 
release of all discussed supernova simulations 



amounts of computer time. First attempts in this direc- 



tion are presently under way (cf. Hoffich & Stein 2002). 

An improvement of the combustion model already 
mentioned in Reinecke et al. (1999b| ), i.e. the full recon- 
struction of all thermodynamic quantities from their mean 
values in every grid cell cut by the burning front, has now 
been compl eted and wa s applied to the Landau-Darrieus 
instability ( Ropke 2002 ). This new model should be im- 
plemented into the full code which, in principle, seems to 
be possible but difficult. In any case, this would increase 
the predictive power and reliability of the models. 

In a similar direction, it might be worthwhile to test 
other subgrid models in combination with our level set 
method. We do not expect major changes because if an- 
other subgrid-model would give us higher (or lower) burn- 
ing speed on the grid scale this would lead to more (or 
less) damping of small-scale structures leaving the prod- 
uct of the two (and therefore the fuel consumption rate) 
approximately unchanged. However, again, such studies 
are presently in progress. 

To conclude, we think that we have made another step 
towards the understanding of type la supernovae. 
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